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PREFACE 


This  document  is  the  final  research  report  on  the  investigation  of  a 
mathematical  algorithm  to  do  target  tracking,  using  doppler  compensated 
correlation  techniques  on  input  time  series  streams  from  several  passive 
acoustic  sensors.  The  algorithm  was  developed  and  programmed  into  a 
testbed  on  a  VAX-750  computer  and  was  tested  using  simulated  time 
series  data  generated  by  the  Tetra  Tech  Broadband  Signal  Simulator. 
Algorithm  performance  proved  dissapointing  due  to:  (1)  Numerical 
instabilities  induced  by  structural  anomolies  in  the  sample  signal 
autocorrelation  function;  (2)  The  extreme  sensitivity  of  objective  function 
to  choice  of  signal  characteristics  and  processing  parameters;  (3) 
Computational  intensity  of  the  algorithm. 
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1.  INTRODUCTION 


For  the  past  year,  Tetra  Tech  has  been  involved  in  the 
development  and  analysis  of  an  algorithm  for  tracking  maneuvering 
submarines  using  Doppler  compensated  correlation  techniques.  The 
intended  goal  was  the  identification  of  an  algorithm  which  would  be 
more  responsive  to  target  kinematic  changes  than  the  usual  Kalman 
filter,  thereby  providing  timely  and  accurate  estimates  of  target 
position  course  and  speed  at  various  times  along  the  track.  The  work 
was  carried  out  under  contract  N00014-8I-C-0535  with  the  Office 
of  Naval  Research. 

The  algorithm  was  developed  and  a  testbed  computer  code  was 
generated  for  implementing  and  testing  the  algorithm.  This  program 
was  written  in  FORTRAN-77  on  a  VAX-750  computer  and  Is  set  up  to 
use  the  Tetra  Tech  Broadband  Signal  Simulator  (BSS)  outputs  as  as 
its  input  time  series.  No  attempt  was  made  in  this  effort  to  use 
actual  at  sea  data.  The  FORTRAN  listings  of  the  testbed  program  are 
included  in  Appendix  A  of  this  report. 

The  algorithm  research  carried  out  by  Tetra  Tech  assumed  the 
availability  of  base  banded,  bandlimited  digitally  sampled  time  series 
from  several  passive  acoustic  sensors.  The  availability  of 
information  such  as  bearings  and  time  delays  were  not  assumed  in 
the  testbed  program,  since  it  was  felt  that  working  directly  with  the 
time  series  data  would  provide  a  more  convenient  means  of 
inplementing  a  maneuvering  target  tracker.  Due  to  the  general 
structure  of  the  algorithm,  additional  measurement  types  such  as 
those  mentioned  above  can  be  easily  added  if  desired. 

As  is  well  known,  one  of  the  characteristics  of  a  sequential  or 
Kalman  type  of  estimating  scheme  is  the  tendency  of  the  estimator 
to  build  up  "inertia"  and  thereby  make  it  unresponsive  to  changes  in 
target  course  and  speed  after  long  periods  of  tracking  time.  This  may 
be  overcome  by  certain  ad  hoc  schemes  such  as  frequent 
reinitializations  or  possibly  by  tampering  with  the  weights  so  as  to 
cause  the  algorithm  to  have  a  shorter  memory,  etc. 


In  view  of  this  it  seemed  reasonable  to  employ  an  estimation 
scheme  which  works  directly  with  selected  blocks  of  time  series 
data  in  a  batch  mode,  and  which  can  be  highly  overlapped  from 
estimation  to  estimation.  Once  the  algorithm  has  been  initialized  and 
is  operation,  the  previous  estimate  of  target  state  can  be  used  to 
initialize  the  trial  solution  for  the  current  estimation.  The 
covariance  matrix  of  the  initialization  state  is  not  carried  over  and 
therfore  the  process  is  without  a  memory.  However,  if  sufficient 
overlap  in  the  input  time  series  is  used,  the  output  states  should 
show  mimfmal  change  from  estimation  to  estimation  while  still 
reflecting  the  most  current  information  available  from  the  time 
series. 

Time  series  generated  by  a  single  moving  source  and  received 
at  two  or  more  spacially  separated  sensors  will  exhibit  different 
Doppler  and  time  delay  characteristics  at  each  of  the  receivers. 
These  characteristics  are,  of  course,  dependent  on  sensor-target 
geometry  and  kinematics,  as  well  as  sound  propagation  physics.  By 
applying  the  appropriate  time  and  Doppler  compensation  to  the 
received  time  series,  pairwise  time  series  correlations  between 
sensors  can  be  maximized.  By  linking  the  time  and  Doppler 
compensation  to  assumed  target  motion,  one  can  adjust  the  target 
state  parameters  to  maximize  (or  minimize)  an  appropriately  chosen 
function  of  the  corresponding  pairwise  cross  correlation  estimates. 

The  assumption  of  digitally  sampled  time  series  sampled  at  a 
uniform  sample  rate  suggests  that  the  time  and  Doppler 
compensation  be  done  In  the  time  domain  using  an  interpolate 
resampling  technique.  This  invloves  estimating  time  series  values 
whose  sample  times  lie  between  the  discreet  sample  times  of  the 
input  time  series.  For  band  limited  signals,  the  well  known  Sampling 
Theorem  provides  a  rational  means  of  performing  the  required 
interpolation  using  neighboring  time  series  points  and  the  sine 
function  as  an  interpolating  function.  This,  in  effect,  generates  a 
piecewise  continuous  representation  of  the  original  time  series 
thereby  allowing  resampling  at  arbitrary  times  which  are  in  concert 


with  trial  target  state  parameters.  Also,  such  a  scheme  allows 
analytic  evaluation  of  the  gradient  vector  with  respect  to  the  state 
variables  and  weighting  vectors. 

As  has  been  mentioned  above,  Tetra  tech  has  implemented 
these  ideas  into  a  testbed  algorithm  on  the  VAX-750  digital 
computer.  Inputs  to  the  algorithm  consist  of  up  to  10  channels  of 
time  series  data.  For  each  channel,  the  algorithm  requires  an 
estimation  of  signal  to  noise  ratio  (SNR)  along  with  the  standard 
deviation  of  the  SNR  for  that  channel.  The  algorithm  is  also  sensitive 
to  such  inputs  as  Integration  time,  processing  bandwidth  and  center 
frequency,  station  location,  sound  speed  in  water,  time  series 
overlap,  and  initialization  parameters  such  as  position  course  and 
speed.  The  target  kinematic  model  consists  of  polynomials  in  water 
time  of  up  to  degree  5  for  x,y  and  z.  The  order  of  the  polynomials  is 
user  specified.  Model  output  consists  of  estimated  target 
parameters,  their  associated  error  variances,  and  the  size  and 
orientation  of  the  2-0  containment  ellipsoid. 

The  algorithm  has  been  tested  using  simulated  time  series 
generated  by  the  BSS.  The  BSS  can  emulate  the  complex  time  series 
generated  by  a  moving  target  having  user  specified  kinematic  and 
spectral  output  signal  characteristics. 

Section  2  of  this  report  gives  an  overall  description  of  the 
workings  of  the  algorithm.  Section  3  describes  the  Gauss-Newton 
estimation  scheme  as  it  applies  to  this  effort,  and  Section  4 
outlines  and  justifies  the  doppler  compensation  scheme  that  we 
have  employed.  Section  5  contains  the  objective  function 
minimization  process  description  and  algorithm  performance  is 
reported  in  Section  6.  Finally,  Section  7  presents  the  summary  and 
conclusions  of  this  effort.  Appendix  A  contains  the  FORTRAN  listings 
of  the  testbed  program  developed  as  part  of  this  effort. 


2. 


ALGORITHM  DESCRIPTION 


The  algorithm  estimation  scheme  assumes  the  availability  of 
several  channels  of  bandlimited,  basebanded,  discreetly  sampled 
digital  data  for  which  all  of  the  pertinent  parameters  such  as  sample 
rate,  bandwidth,  and  center  frequency  are  known,  and  all  of  which 
contain  signal  from  a  common  emitter.  For  the  purposes  of  this 
analysis,  we  will  assume  that  the  noise  on  each  channel  is  mutually 
uncorrelated.  We  will  also  assume  that  the  target  is  moving  along 
some  3-dlmensional  trajectory,  which  is  given  by  the  vector  function 
P(s;t),  where  s  is  the  state  vector  to  be  estimated  and  t  is  the  time 
along  the  trajectory. 

The  testbed  version  of  the  algorithm  assumes  that  each  of  the 
components  of  P(s;t)  is  an  n'th  order  polynomial  in  t,  and  the  state 
vector  s  consists  of  the  coefficients  of  these  polynomials.  The 
testbed  user  may  specify  n  to  be  any  non  negative  integer  up  to  and 
including  5.  In  practice,  n  is  usually  chosen  to  be  1,  resulting  in 
linear  target  motion  at  constant  speed.  In  this  case,  P(s;t)  is  given 
by 

P(s;t)=Po+Vt  (2.1) 

and  the  state  vector  s  may  be  represented  in  transposed  form  by 

s=[P0TVTlT  (2.2) 


The  superscript  T  denotes  the  usual  matrix  transpose  operator. 

In  order  to  keep  the  algorithm  tractible,  we  have  assumed 
linear,  constant  speed  sound  propagation.  More  complicated 
propagation  models  could  have  been  incorporated,  but  it  was  deemed 
an  unnecessary  complication  at  these  early  stages  of  algorithm 
development  and  feasibility  analysis. 

The  central  idea  of  the  algorithm  is  to  mimimize  a  quadratic 
form  of  "system  functions".  The  system  functions  are  dependent  on 
the  collection  of  pairwise  normalized  sample  correlation  envelope 
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functions  which  have  been  adjusted  to  account  for  assumed  target 
kinematics.  Each  sample  correlation  envelope  function  is  obtained  by 
correlating  the  samples  from  a  selected  reference  channel  with 
modified  sets  of  samples  that  have  been  interpolated  and  resampled 
from  each  the  other  channels  comprising  the  tracking  system.  The 
resampling  times  are  calculated  as  a  function  of  the  current  value  of 
the  state  vector,  sensor  kinematics,  the  assumed  propagation  model, 
and  channel  signal  processing  parameters  such  as  center  frequency, 
sample  rate,  and  bandwidth. 

To  be  specific,  let  us  consider  a  pair  of  channels,  say  channels 
X  and  Y.  Pick  a  set  of  samples  {X1,X2,...,Xn}  from  channel  X.  These 
samples  correspond  to  arrival  times  {ut,u2,...,un}  on  channel  X.  In 
order  to  do  the  motion  compensation,  we  need  to  calculate  the 
corresponding  arrival  times  {v1,v2,...,vn}  on  channel  Y.  This  is  done  by 
using  the  candidate  source  trajectory  P(s;t)  to  calculate  emitter 
times  {t|,t2,...,tn}  corresponding  to  the  X  channel  arrival  times 
{ui,U2,...lun}l  and  using  these  emitter  times  to  project  the 
corresponding  arrival  times  {v1,v2,...,vn}  on  channel  Y.  We  then 
interpolate  and  resample  the  Y  channel  at  the  {vt,v2,...,vn}  thereby 
obtaining  a  new  set  of  samples  {Yj,Y2,...,Yn}.  The  interpolation  is 
accomplished  using  a  truncated  sine  function  as  an  interpolating 
function.  Details  of  the  interpolation  scheme  are  provided  in  Section 
4  of  this  report. 

The  magnitude  squared  cross  correlation  estimate  tfxy  is  then 
calculated  by 

«xH2xivi*IA2l><il2Zlvil2}  (2-3) 

where  the  three  sums  in  the  above  expression  are  taken  over 
1=1,2,...,n  and  the  superscript  (*)  denotes  complex  conjugation.  The 
values  of  the  tfxy  so  obtained  are  used  to  form  the  aformentioned 
system  functions  Fxy  which  are  used  in  the  minimization  process. 
The  system  functions  are  given  by 

Fxy=ln(Gxy/tfxy)  (2.4) 
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where  Gxy  is  the  a  priori  expected  value  of  tfxy.  Gxy  Is  related  to  the 
Input  SNR  estimates  on  each  channel  by 

Gxy=((  1  +SNRx~1)(  1  +SNRy-,)]_1  (2.6) 

Note  that  the  Fxy  are  chosen  such  that  they  have  value  0  when  given 
error  free  information. 

Finally,  Q(s),  a  positive  definite  quadratic  form  of  the  Fxy,  is 
formed  ov'T  all  channel  pairs,  and  by  using  gradient  methods,  the 
state  vector  s  is  adjusted  so  as  to  mimimize  Q(s).  The  vector  s0 
which  mimimizes  Q(s)  is  taken  as  the  state  estimate. 

A  fallout  of  the  minimization  process  is  an  estimate  of  the 
state  covariance  matrix.  This  matrix  is  used  to  calculate  the 
ellipsoidal  containment  region  which  provides  the  user  with  a 
geometric  indication  of  algorithm  performance. 


3.  THE  ESTIMATION  SCHEME 


3.1  Preliminary  Details  and  Notation 

this  section  details  the  estimation  scheme  in  rather  general 
mathematical  terms.  Preliminary  to  the  discussion  we  establish  the 
following  notation  which  will  be  used  throughout  the  remainder  of 
this  report. 

If  X  and  Y  are  n-dimensional  complex  valued  vectors,  the 
complex  inner  product  of  X  and  Y,  denoted  by  <X,Y>,  is  defined  by 

<X,Y>=^X1Yj*  (3.1) 

where  the  sum  is  taken  over  i=1,2,...,n,  and  the  X^  and  Yj  are  the 
complex  valued  components  of  X  and  V,  respectively.  The  (*) 
notation  denotes  complex  conjugation.  Note  that  the  complex  inner 
product  is  conjugate  symmetric  in  that  the  following  relationship 
holds: 


<Y,X>=<X,Y>* 


(3.2) 


We  define  the  norm  of  X,  denoted  by  |X|,  by 


|X|=/<X,X> 


(3.3) 


In  this  notation,  Equation  2.3  of  Section  2  becomes 

tfxy  =  |  <X,Y>  1 2/{  |Xj|2  j|Y|2}  (3.4) 

If  each  of  the  components  of  the  complex  valued  vector  is  a 
differentiable  function  of  some  real  parameter  0,  then  denote  the 
vector  consisting  of  the  corresponding  derivatives  (partial 
derivatives)  by  dX/d8  (9X/08).  It  is  easy  to  verify  that  the  following 
useful  relationship  is  true 
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The  following  figures  show  the  kinds  of  difficulties  one  faces 
even  under  the  best  of  circumstances.  Figures  6.2-6.5  show  plots  of 
the  objective  function  versus  errors  in  the  probe  or  trial  solution 
velocities(ft/sec)  in  both  the  x  and  y  direction  for  several  sets  of 
processing  parameters.  The  assumed  position  error  is  assumed  to  be 
fixed  at  0  ft.  The  velocity  errors  range  from  ±3  ft/sec  in  both  the  x 
and  y  direction.  The  center  of  the  portion  of  the  x-y  plane  shown 
represents  0  error.  The  2  axis  represents  the  valuer  of  the  objective 
function  and  has  been  inverted  for  the  sake  of  this  presentation.  The 
z  axis  ranges  from  0  (top)  to  1600  (bottom).  The  signal  bandwidth  is 
assumed  to  be  8  Hz. 

Figure  6.2  was  generated  assuming  a  center  frequency  of  20  Hz 
with  an  integration  time  of  20  seconds  and  presents  a  very  clean 
objective  function  over  the  plot  range.  In  Figure  6.3,  the  integration 
time  has  been  increased  100  seconds.  The  resulting  plot  exhibits  a 
much  more  spiked  peak  around  (0,0)  and  the  outskirts  show  a  good 
deal  of  ripple.  Figure  6.4  is  similar  to  Figure  6.2  except  that  the 
center  frequency  of  the  signal  has  been  shifted  to  250  Hz.  This  band 
shift  has  caused  the  peak  to  become  very  sharp  with  even  more  ripple 
evidenced  in  the  outskirts.  The  processing  situations  depicted  in 
Figures  6.3  and  6.4  would  present  problems  to  the  tracker.  Figure  6.5 
has  identical  parameters  as  Figure  6.2  with  the  exception  that  the 
assumed  position  error  has  been  offset  to  300  ft.  in  both  the  x  and  y 
directions  resulting  in  a  much  more  planar  shape  and  having 
increased  the  overall  magnitude  of  the  objective  function 
considerably. 

Figures  6.6  and  6.7  show  plots  of  the  objective  function 
versus  errors  in  the  trial  position  estimates  with  the  geometry  in 
Figure  6.1  applying.  The  position  errors  range  from  +  400  ft  in  both 
the  the  x  and  y  directions.  Figure  6.6  corresponds  to  a  signal 
bandwidth  of  8  Hz  and  presents  a  rather  smooth  function  with 
essentially  no  unusual  structure.  In  Figure  6.7  we  have  increased  the 
signal  bandwidth  to  100  Hz  and  moved  the  center  frequency  to  80  Hz. 
The  resulting  plot  contains  considerable  structure  with  the 


parameters,  we  developed  a  computer  program  to  generate  the 
expected  value  of  Q(s)  for  a  particular  signal  autocorrelation 
function.  The  outputs  are  sensitive  to  the  aformentioned  signal  and 
processing  parameters  as  well  as  target/sensor  geometry  and 
kinematics.  The  autocorrelation  function  we  have  chosen  is 
triangular  on  the  interval  I-T,T1  and  has  a  spectral  density  function 
of  the  form  sinc2(o>T/2).  Its  bandwidth  is  approximately  1/T.  This 
particular  autocorrelation  function  has  the  advantage  of  being 
integrable  in  closed  form.  The  program  is  setup  to  generate  surfaces 
of  expected  values  of  the  tracker  objective  function,  holding  the 
position  probes  fixed  and  letting  the  velocity  probes  vary,  or  holding 
the  velocity  probes  fixed  and  letting  the  position  probes  vary. 

Some  results  are  presented  for  the  geometry  shown  in  Figure 
6.1.  Here  the  target  is  assumed  to  be  in  the  center  of  a  square  box 
with  the  sensors  located  at  the  vertices.  The  distance  from  the 
target  to  each  of  the  sensors  is  assumed  to  be  2  nmi.  A  SNR  of  6  dB 
with  a  <*5nr=3  dB  is  assumed. 


o 
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Figure  6. 1 

Test  Geometry 
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ALGORITHM  PERFORMANCE 


6.1  Parameter  Selection 

Since  the  objective  function  Q(s)  is  a  very  complicated 
function  (most  of  which  is  carried  out  in  the  complex  domain)  of  a 
number  of  processing  parameters,  some  thought  had  to  be  given  to 
their  effects  on  algorithm  behavior.  The  user  has  control  over  such 
things  as  processing  bandwidth,  center  frequency,  and  integration 
time. 

For  example,  if  one  chose  to  process  a  signal  having  a 
sufficiently  high  center  frequency  for  a  sufficiently  long  period  of 
time,  Q(s)  could  become  quite  sensitive  to  trial  solution  errors  in 
velocity.  In  fact  one  would  expect  to  see  spike  like  behaviour  in  the 
objective  function  in  the  neighborhood  of  the  true  velocity 
components,  which  could  conceivably  cause  convergence  difficulties. 

Similarly,  a  wide  bandwidth  signal  could  cause  similar  spike 
like  behavior  in  trial  solution  position  estimates.  Small  positional 
errors  wind  up  in  the  ripple  of  the  objective  function  which  wreaks 
havoc  on  the  convergence  process  we  have  chosen. 

On  the  other  hand,  if  the  processing  band  is  too  narrow,  the 
lack  of  time  delay  resolution  may  not  provide  any  meaningful 
information  to  the  algorithm,  again  resulting  in  poor  behavior. 

These  were  in  fact  some  of  the  problems  that  were 
encountered  during  the  early  stages  of  algorithm  testing.  Since  the 
testbed  computer  code  was  newly  developed,  it  was  not  known 
whether  poor  early  algorithm  performance  was  due  to  bugs  in  the 
program,  or  whether  the  algorithm  just  did  not  work,  or  whether  we 
had  just  chosen  bogus  processing  parameters.  After  considerable 
reexamination  , rechecking,  and  rederiving  the  mathematics,  we 
decided  that  they  were  correct.  We  could  find  no  bugs  in  the  program 
and  so  our  only  recourse  was  to  give  careful  scrutiny  to  our  choice 
of  test  parameters. 

In  order  to  gain  insight  to  the  effects  of  processing 


which  yields 


ds=-(FsT WFS)’ 1  (Fst WF)=Css  (Fst WF)  (5.5) 

The  advantage  of  this  approach  is  that  it  does  not  require  the 
computation  of  second  derivatives,  and  that  an  estimate  of  the  state 
covariance  matrix  falls  out  of  the  process. 

In  summary,  if  s  is  the  current  trial  solution  for  the  process, 
we  form  the  next  iterate  s'  by  calculating  9s  as  in  the  preceeding 
equation  and  forming  s'  by 


s'=s  +  9s  (5.6) 

To  stop  the  process  we  check  to  see  if  the  magnitude  of  9s  is  within 
some  predescribed  tolerance.  If  so  the  process  is  stopped  and  the 
current  value  of  the  state  vector  is  returned  as  the  solution.  If  not 
the  process  continues  until  a  solution  is  returned  or  the  maximum 
iteration  count  is  exceeded. 

The  entire  process  is  described  in  Figure  5.1. 


5.  THE  MINIMIZATION  PROCESS 


5.1  The  Iteration  Scheme 

The  function  Q(s)  is  a  complicated  function  of  the  state  vector 
and  measurement  vector,  the  mimimization  of  which  does  not  seem 
amenable  to  closed  form  solutions.  We  therefore  must  rely  on 
iterative  techniques  to  solve  the  problem.  The  following  paragraphs 
describe  the  technique  we  have  chosen  to  accomplish  the 
minimization. 

Recall  that  Q(s)  is  a  quadratic  form  in  the  system  functions 
and  may  be  written 

Q(s)=FtWF  (5,1) 

where  the  weighting  matrix  W  is  chosen  to  be  the  inverse  of  the 

covariance  matrix  of  the  system  residual  vector.  Under  the 
assumption  of  slowly  varying  weights,  the  gradient  vector  of  Q(s) 
may  be  written 

VQ(s)=2FstWF.  (5.2) 

At  a  local  mimimun  we  have  the  necessary  condition  that 

VQ(s)=2FstWF=0  (5.3) 

An  approach  which  has  been  used  sucessfully  is  demonstrated  in  the 

following  discussion.  Suppose  that  the  algorithm  has  reached  a  stage 
such  that  s  is  the  current  value  of  the  trial  state  vector.  We  would 
like  to  find  the  perturbation  3s  to  add  to  s  which  will  improve  the 
estimate.  A  reasonable  approach  is  to  solve  the  the  following 
perturbed  gradient  equation  for  3s 


F8tW(F  ♦  F83s)=0 


(5.4) 


receiver  times  using  Q^(t)  as  the  interpolating  function,  thereby 
generating  a  set  of  Y  channel  samples  which  reflect  the  Doppler 
corrections  Implied  by  the  current  value  of  the  state  vector. 

In  order  to  obtain  the  receiver  times  on  the  Y  channel,  we 
solve  the  following  pair  of  equations  for  vk  given  uj<: 

uk=Tk  ♦  |  P(s,Tk)-px  |  /c  (4.6) 

W  j  P(s,Tk)-Py  j  /c 

where  Px  and  Py  are  the  respective  position  vectors  of  the  X  and  Y 

channel  receivers,  c  is  the  speed  of  sound  in  water,  and  is  the 
emitter  time.  This  is  done  by  solving  the  first  equation  for  Tj<  using  a 
Newton-  Raphson  technique,  and  then  using  the  second  equation  with 
the  value  of  so  obtained  to  obtain  v^. 

Recall  that  we  are  assuming  that  all  of  the  input  time  series 
data  has  been  basebanded  from  some  center  frequency  fc.  Because  of 
this,  a  phase  correction  prior  to  correlation  is  required  on  both 
channels.  This  merely  amounts  to  heterodyning  the  samples  back  up 
to  their  original  center  frequency  fc#  Form  the  complex  vectors  X  and 
Y  whose  k'th  components  are  given  by  exp{2Ttjfc(uk-u0)  >Xn+k_ !  and 
exp{2TT jfc(vk-v0)  }Yn+|<_j,  respectively.  Then  tfxy  is  given  by 

tfxy=  |  <X,Y>  1 2/{|X|2I  V|2>  (4.7) 

4.3  System  Function  Derivatives 

The  algorithm  uses  gradient  methods  to  minimize  Q(s),  which 
is  the  quadratic  form  of  the  system  functions  Fxy  as  discussed 
above.  This  necessitates  the  calculation  of  the  derivatives  of  Q(s) 
with  respect  to  each  of  the  state  variables.  Most  of  the  work  is  done 
in  computing  the  derivatives  of  the  tfxy  with  respect  to  the 
measurement  vector  and  each  of  the  state  variables.  The 
mathematical  development  of  these  derivatives  is  straightforward 
but  extremely  tedious  and  will  not  be  included  here. 


reconstructed  from  its  samples,  provided  that  the  digital  sample  rate 
is  at  least  as  great  as  the  bandwidth  of  the  signal.  The 
reconstruction  uses  the  "sine  x  over  x"  or  sine  as  an  interpolating 
function.  Specifically,  if  Z(t)  is  a  time  series  having  non  zero 
frequency  content  only  in  the  interval  [-b/2,b/2l,  and  if  Z(t)  is 
uniformly  sampled  over  all  time  at  a  sample  rate  f(j>b,  producing 
samples  Zn,  -oosnsoo,  and  such  that  Z0  corresponds  to  a  sample  time 
of  0,  then  Z(t)  is  reconstructed  exactly  from  its  samples  by 

Z(t)=£sinc{TT(fdt-n)}Zn  (4.3) 

The  above  sum  is  taken  over  all  n.  This,  however,  involves  summing 
over  an  infinite  number  of  elements.  It  therefore  seems  reassonable 
to  approximate  the  reconstruction  of  Z(t)  from  its  samples  by  using 
a  time  limited  version  of  the  sine  function.  If  we  define  the 
interpolating  function  Qh(t)  by 

sinc(t),  1 1 1  sh 

Qh(t)=  (4.4) 

0  | t | >h 

,  then  Z(t)  may  be  approximated  by 

Z<t)*2QhMfdt"n)>zn  (4.5) 

which,  for  any  given  value  of  n  involves  only  finite  sums. 

Let  us  assume  we  are  working  with  channels  X  and  V  and  we 
wish  to  calculate  tfxy  for  the  current  value  of  the  state  vector.  Pick 
a  set  of  reference  samples  {Xn,Xn+1,...,Xn+^-1}  from  channel  X.  These 
samples  correspond  to  the  set  of  channel  X  receiver  times 
{un,un+1,...,un+tt-ih  We  the  use  P(s;t)  and  the  assumed  propagation 
model  to  determine  the  corresponding  set  of  channel  V  receiver 
times  {vn,vn+),...,vn+^j_1}.  Channel  Y  is  then  interpolated  at  these 


DOPPLER  COMPENSATION  SCHEME 


The  algorithm  assumes  the  availability  of  M  channels  of  data, 
and  that  the  time  series  for  each  data  channel  contains  basebanded, 
bandlimited  data  that  has  been  sampled  at  at  least  the  Nyquist 
sampling  rate.  For  the  sake  of  notational  simplicity,  we  will  assume 
that  all  of  the  channels  have  the  same  center  frequency,  fc,  and  the 
same  digital  sampling  rate  fd.  We  also  assume  that  the  n'th  sample 
on  each  channel  occurs  at  the  same  time.  The  time  between  samples, 
At,  is  given  by 


At=1/fd 


(4.1) 


Therefore  if  the  O'th  sample  corresponds  to  t0,  then  the  n'th  sample 
corresponds  to  tn,  where 


tn=t0  +  nAt 


(4.2) 


The  algorithm  requires  an  estimate  of  the  signal  to  noise  ratio 
and  its  corresponding  error  variance  on  each  channel  of  input  data.  It 
is  assumed  that  these  are  independently  specified  and  will  be  denoted 
SNRm,  where  the  subscript  m  refers  to  the  particular  channel 
designator. 


In  order  to  effect  the  proper  Doppler  compensation,  it  is 
necessary  to  interpolate  between  samples  in  the  time  domain,  with 
the  interpolation  times  reflecting  the  current  value  of  the  state 
vector. 

The  well  known  sampling  theorem  from  signal  processing 
states  that  a  complex  valued  bandlimited  signal  can  be  completely 


Suppose  we  want  to  add  a  new  measurement  set  to  the 
algorithm  that  is  independent  of  those  already  incorporated  and  can 
be  described  with  a  single  system  equation.  This  merely  requires  the 
specification  of  the  corresponding  system  function  and  its 
associated  partial  derivatives.  In  keeping  with  our  earlier 
discussion,  let  us  further  assume  that  the  partial  of  the  new  system 
equation  with  respect  to  its  measument  vector  are  independent  of 
the  state  vector.  In  this  case  we  first  form  the  new  scalar  weight  w 
by 

w=l/3(aF/dmi)Omi)2  (3.19) 

where  F  is  the  new  system  function,  the  m^  are  its  associated 
measurements  and 

'  dF/arn^IdF/am*].  (3.20) 

The  updated  B  vector  and  Z  matrix  are  given  by 

B=Bold  +  wF(9F/aS)  (3.21) 

z=z0id  +  w  (aF/as)  (aF/as)T  (3.22) 

where  3F/a5  is  the  vector  of  partlals  of  F  with  respect  to  the 
elements  of  the  state  vector. 


where  Fs  denotes  the  matrix  whose  mn'th  element  is  3Fm/3sn>  If  the 
above  equation  holds,  a  perturbation  in  the  measurement  vector  3m 
induces  a  perturbation  in  the  state  vector  3s,  which  to  within  first 
order  terms,  obeys  the  relationship 

FsT  W(F  +  Fs3s  +  Fm3m)=0.  (3. 1 4) 

This  implies 

3s*-(Fst  WFs)_1Fst  WFm3m,  (3. 1 5) 

which  yields,  after  substituting  Cmm_1  for  W 

CSs%^^s3st)=(FsTWFs)"1  (3.16) 

The  expression  for  Css  given  above  is  extremely  convenient  and  can 
be  used  to  provide  geometrical  insight  to  algorithm  behaviour.  In 
particular,  it  is  used  to  derive  the  ellipsiodal  containment  region  of 
the  current  estimate  of  the  state  vector. 

3.3  Incorporation  of  Additional  Measurements 

The  form  of  the  estimator  used  in  this  algorithm  has  the 
advantage  that  it  Is  easy  to  Incorporate  new  types  of  measurements 
should  the  need  arise.  If,  for  example,  the  system  can  provide 
independent  estimates  of  position  and  velocity  or  if  another 
independent  sensor  comes  on  line,  the  new  measurements  so  provided 
may  be  entered  as  an  additive  partitions  to  the  system  weighting 
matrix  and  gradient  vector.  Let  us  first  establish  the  following 
notation.  Let 


Z=  (FstWF«) 


(3.1?) 


We  now  provide  our  rational  for  our  choice  of  W.  Let  m  denote 
the  vector  of  measurements  associated  with  F.  In  our  case,  each 
component  of  m  is  an  SNR  estimate  from  one  of  the  output  channels. 
For  any  fixed  value  of  s,  the  perturbation  in  the  vector  F  induced  by  a 
perturbation  In  the  vector  m  is  approx  .rated  by 

^FMi-’mldm  (3.9) 

Here  lFm]  is  the  matrix  whose  Ij'th  element  is  dF^Smj.  Assuming  the 
error  distribution  is  zero  mean  with  covariance  matrix  Cmm,  and 
that  the  approximation  to  0F  holds  over  the  probable  values  of  m,  we 
may  approximate  the  covariance  matrix  of  F  by 

^FF**rmcmm,:mT  (3.10) 

The  weighting  matrix  W  referred  to  in  in  Equation  (3.8)  is 
chosen  to  be  (Cpp)-1.  Therefore  the  scalar  function  Q(s)  is  given  by 

Q(s)=FT(CFFr'F  (3.11) 

This  choice  of  W  has  intuitive  appeal  in  that  measurements  with 
higher  variance  get  less  weight  and,  therefore,  have  less  of  an  effect 
on  the  final  outcome. 

We  now  turn  our  attention  to  the  approximation  of  the  output 
state  covariance  matrix,  which  will  be  denoted  by  CgS.  The 

components  of  the  state  vector  s0  that  minimizes  Q(s)  satisfy  the 
equation 

dQ/0Sj=2(dF/dSj)WF  +  F^dW/as^O  ,i=1,2 . k  (3.12) 

If  we  assume  slowly  varying  weights  which  enables  us  to  ignore  the 
second  term  in  the  above  equation,  then  s0  satisfies  the  system  of 
equations 


d<x,Y>/ae=<ax/ae,Y>  ♦  <x,aY/ae>.  (3.5) 

The  above  relationship  is  indespensible  in  computing  the  gradient 
derivatives  during  the  minimization  process. 

3.2  Theoretical  Development 

The  general  theoretical  basis  for  the  estimation  scheme  is  an 
extension  of  the  techniques  of  linear  optimal  estimation  theory  to 
the  nonlinear  case.  We  evaluate  a  set  of  system  equations  Fxy  as 
given  by  Equation  ( )  and  then  minimize  a  positive  definite  quadratic 
form  of  the  Fxy.  Suppose  we  have  k  such  system  functions  for  a 
particular  application.  Changing  notation  for  convenience,  let  the 
system  functions  be  denoted  by  F1,F2,...,F((and  let 

F(s)=[F1,F2,...,  Fk]T  (3.6) 

denote  the  k-dimensional  vector  of  system  equations.  Note  that  the 
solvability  of  the  above  equations  implies  that  ngSk,  where  ns  is  the 
dimentionality  of  the  state  vector. 

If  the  geometric  and  signal  assumptions  are  perfectly 
compatible  with  the  measurement  data,  there  exists  a  value  of  the 
state  vector  s0  for  which 

F(s)=0.  (3.7) 

In  general,  however  the  data  does  not  support  perfect 
compatibility,  and  for  each  value  of  the  state  vector  s,  F(s)  may  be 
considered  a  vector  of  F-residuals.  An  optimal  estimate  of  s  is  a 
vector  which  minimizes  a  particular  quadratic  form  in  the 
F-residuals.  If  W  is  an  appropriately  chosen  positive  definite 
symmetric  matrix,  the  optimal  estimate  is  the  value  which 
minimizes  the  scalar  function 


Q(s)=FtWF. 


(3.8) 


crisscrosses  being  corresponding  to  time  delays  among  the  various 
sensors.  Again,  finding  the  minimum  of  such  a  function  if  the 
solution  starts  off  of  the  main  peak  presents  quite  a  formadlble  task 
to  the  tracker. 

We  have  presented  these  results  to  demonstrate  difficulties  in 
choosing  a  set  of  operating  parameters  for  which  we  may  have  some 
hope  of  achieving  positive  results.  To  this  end  and  after  having  a 
large  number  of  cases  as  above,  we  chose  to  work  with  a  time  series 
centered  at  20  Hz  containing  an  8  Hz  signal  in  10  Hz  wide  total 
processing  bin.  The  signals  were  generated  by  the  BBS  for  sensors 
having  the  geometry  described  in  Figure  6.1  above,  and  had  an  overall 
SNR  of  +6dB  in  the  processing  band  zo  each  of  the  sensors.  Using  an 
integration  time  of  20  seconds  and  initialiazing  the  algorithm  with 
"truth"  at  various  points  along  the  target  track,  the  algorithm 
evidenced  unstable  convergence  in  nearly  every  case.  Since  we  had 
gone  over  the  computer  program  and  the  mathematics  very  carefully 
and  could  find  no  errors,  we  had  to  look  elsewhere  for  an  explanation 
of  the  poor  algorithm  performance.  We  believe  the  answer  lies  in  the 
form  of  the  system  functions. 

The  Z  matrix  and  B  vector  defined  in  Section  3  of  this  report 
are  used  extensively  in  the  Gauss-Newton  mimimization  technique 
that  we  chose  to  implement  for  this  algorithm.  Note  that  both  Z  and 
B  contain  the  partial  derivatives  of  the  system  functions  with 
respect  to  each  of  the  elements  of  the  state  vector.  When  the  system 
is  near  a  solution,  the  partials  of  the  system  functions  are  near  0, 
since  they  are  essentially  the  derivatives  of  the  signal 
autocorrelation  function  at  t=0.  This  causes  the  Z  matrix  to  become 
numerically  unstable  as  the  process  nears  a  solution,  thereby  causing 
unpredictable  algorithm  behaviour.  Efforts  to  accomodate  alternate 
forms  of  the  system  functions  turned  proved  unsucessful  because  we 
could  not  find  a  tractible  method  of  comparing  time  series  that  did 
not  involve  correlations.  The  only  other  alternative  to  sucessfully 
implement  the  algorithm  would  be  to  incorporate  a  more 
sophisticated  minimization  technique  which  would  overcome  the 


instability  problems.  This,  too,  could  have  drawbacks  since  most 
such  algorithms  involve  several  one  dimensional  searches  which 
require  several  objective  function  evaluations  which  is  quite 
computationally  intensive. 

We  actually  tried  a  method  which  involves  halving  the  length 
of  the  search  vector  9s  until  Qj<+t<Qk  at  which  time  we  would 
recalculate  a  new  search  vector  and  continue  the  process.  This 
proved  to  be  too  computationally  intensive  even  for  TW  products  on 
the  order  of  200  and  including  4  sensor  geometry.  We  did,  however 
obtain  some  success  in  achieving  algorithm  convergence.  However, 
each  evaluation  of  Q(s)  took  about  2  minutes  CPU  time  on  the  VAX, 
resulting  in  enormous  total  algorithm  processing  times.  This  was 
deemed  unsatisfactory  from  the  point  of  view  of  any  practical 
application. 


7.  SUMMARY  AND  CONCLUSIONS 


An  algorithm  to  use  doppler  compensated  resampling  and 
correlation  techniques  on  digitally  sampled  time  series  was 
developed  and  tested  using  synthetic  time  series  data  generated  by 
the  Tetra  Tech  Broadband  Signal  Simulator  (BBS).  The  algorithm  was 
designed  to  provide  timely  track  estimates  by  using  highly 
overlapped  time  series  segments  from  the  receiving  sensors  in  a 
batch  mode,  thereby  giving  the  process  a  short  memory  and  increased 
responsiveness  to  target  maneuvers.  The  target  motion  model 
consisted  of  polynomial  functions  of  time  of  arbitrary  order  up  to  5. 

Algorithm  performance  proved  dissapointing  for  several 
reasons: 

(1)  Computational  intensity  was  more  than  was  originally 
envisioned. 

(2)  Structure  of  the  correlation  functions  induced  numerical 
instabilities  in  the  convergence  process  which  could  not  be  easily 
overcome. 

(3)  The  structure  of  the  objective  function  is,  in  general,  quite 
irregular,  thereby  requiring  careful,  and  perhaps  limited,  choice  of 
processing  parameters  in  order  to  have  any  hope  of  successful 
performance. 

Some  improvement  in  running  time  could  be  achieved  by 
simplifying  the  resampling  technique  to  use  first  and  perhaps  second 
order  time  series  stretches  based  on  the  current  value  of  the 
position  and  velocity  estimates. 

Curing  the  numerical  instabilities  seems  to  be  the  most 
difficult  hurdle  to  overcome  due  to  the  often  unusual  structures 
evidenced  in  correlation  functions.  For  this  reason,  we  feel  that  any 
further  endea vers  to  improve  upon  such  an  algorithm  would  prove  to 
be  dissapointing  and  recommend  that  further  research  be 
discontinued. 
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